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Abstract 

We consider a particle system on Z d with real state space and interactions of 
infinite range. Assuming that the rate of change is continuous we obtain a Kalikow- 
type decomposition of the infinite range change rates as a mixture of finite range 
change rates. Furthermore, if a high noise condition holds, as an application of this 
decomposition, we design a feasible perfect simulation algorithm to sample from the 
stationary process. Finally, the perfect simulation scheme allows to forge an algorithm 
to obtain an explicit construction of a coupling attaining Ornstein's d— distance for 
two ordered Ising probability measures. 
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1 Introduction 



In this paper we present a Kalikow-type decomposition for interacting multicolor systems 
on Z d having real state space and interactions of infinite range. By a Kalikow-type de- 
composition we mean a representation of the infinite range rates as a countable mixture of 
local change rates of increasing range. This decomposition extends the notion of random 
Markov chains to interacting particle systems and has many potential theoretical conse- 
quences and applications. As a first example we present a perfect simulation algorithm 
which is based on the decomposition. As a corollary we obtain a result about the existence 
and uniqueness of the invariant measure of the system as well as a rate of convergence to 
stationarity. As second application we construct a coupling attaining the d-distance for 
two ordered Ising probability measures. 
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13.081-970 - Campinas, SP - BRAZIL, e-mail addresses: galves@usp.br, nancy@ime.unicamp.br, 
eva.loecherbach@u-cergy.fr and orlandi@mat.uniroma3.it 
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By a perfect simulation algorithm we mean a simulation which samples in a finite window 
precisely from the stationary law of the infinite process. More precisely, for any finite set 
of sites F we want to sample the projection of the stationary law on F. Our approach is 
feasible in the sense that it stops almost surely after a finite number of steps. It does not 
require any duality or monotonicity properties. We do not assume that the system has a 
dual, or is attractive, or monotone in any sense. Our system is not spatially homogeneous. 
The basic assumptions are the continuity of the infinite range change rates together with 
a high- noise condition (Condition f)5.33|) : fast decay of the range influence on the change 
rate and a certain subcriticality-criterion). 

Concerning possible applications, perfect simulation of infinite range continuous (or dis- 
crete) systems has shown to be an important tool, for instance in statistical inference 
for Gibbs distributions, for Bayesian statistics, for maximum likelihood estimation, rates 
of convergence of estimators. This field of research has enormous relevance due to its 
applications in image processing, spatial statistics, gene expression, to cite just a few of 
the possible applications. Some of the applications use discrete state space whereas oth- 
ers deal with continuous state space, some deal with nearest neighbor interaction (image 
processing) whereas others need infinite range interactions (tumor growth, dynamics of 
populations of neurons, or gene expression). We refer the reader to the books of M0ller 
and Waagepetersen (2004), Gaetan and Guyon (2010) for some applications. 

Let us stress that from an applied point of view, it is important to be able to deal with 
infinite range systems. In neuroscience for example, for all practical purposes, the inter- 
actions between neurons have infinite range, see for example Cessac, Rostro, Vasquez and 
Vieville (2009) and Cessac, Nasser and Vasquez (2010). Also, very simple error structures 
superposed on finite range systems will produce infinite range systems. This is the case 
for example of blurred images, see for example Nishimori and Wong (1999) and Tanaka 
(2002). Up to now, denoising of images has been studied mostly with very simple models 
(finite lattice, nearest neighbor interaction) due to the difficulty to sample from a Gibbs 
distribution with infinite range interactions and/or continuous state space. 

So it is important to build up a unifying approach that enables us to deal with infinite 
range interactions and/or continuous state space at the same time. Although perfect 
simulation of Gibbs random fields are known for discrete state spaces and finite range 
interaction, very little is known for simulation of continuous state spaces and/or infinite 
range interaction. To the best of our knowledge ours is the first concrete result concerning 
perfect simulation for interacting particle systems with continuous state space and infinite 
range interactions. 

This paper is organized as follows. The model and the Kalikow-type decomposition (The- 
orem [1]) are presented in Section 2. The aim of Section H] is to present some examples 
where the decomposition can be explicitely done. In particular, we apply Theorem Q] to 
the important case of Gibbs measures with infinite range interactions and continuous spin 
values. In Section [5] we present the perfect simulation algorithm as a main application of 
the Kalikow-type decomposition. In particular, Theorem [3] shows that the proposed algo- 
rithm is feasible under a high noise condition. The proofs are given in Sections [3] and 
In Section \7\ we give as a main application of the perfect simulation algorithm an explicit 
construction of a coupling attaining Ornstein's d— distance for two ordered Ising probabil- 
ity measures. We present a section discussing the user impatient bias. We conclude the 
article with final comments and a bibliographical discussion. 
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2 Definitions, notation and convex decomposition 



We consider interacting particle systems on Z rf having state space A and interactions of 
infinite range. The elements of the state space A are called colors. To each site in Z rf we 
assign a color. The coloring of the sites changes as time goes by. The rate at which the 
color of a fixed site i changes from a color a to a new color b is a function of the entire 
configuration and depends on b. 

In what follows, we suppose that A is a Borel subset of R, equipped with its Borel sigma- 
field A. g will be a finite non-negative reference measure on (A, A). The initial lowercase 
letters a, b, c, . . . will denote elements of A. We denote by S = A z the configuration space 
with its product sigma algebra, <S. We call an element of S a configuration. Configurations 
will be denoted by Greek letters rj, ... A point i G Z rf will be called a site. We define 
on Z d the L 1 norm, ||i|| = Y2k=i ^ or ^ — ^' ^ ^ ne ball °f radius k be 

V t (k) = {j eZ d ;\\j -i\\ <k}. 

As usual, for any i G Z d , r](i) will denote the value of the configuration r\ at site i. 
By extension, for any subset V C Z d , rj{V) G A v will denote the restriction of the 
configuration r] to the set of positions in V. For any r], i and a, we shall denote rf' a the 
modified configuration 

rf' a {j) = ^0)) f° r all i 7^ and rf' a {i) = a. 

For any i G Z d , let Ci(a,r]) be a positive *4<X>tS — £?(IR+) — measurable function such that the 
following two properties hold. Firstly, for £— almost all a G A, r] h-> Cj(a, 77) is continuous. 
Secondly, we have 

sup sup / q(a, r))g(da) < 00. (2-1) 

A multicolor system with interactions of infinite range is a Markov process on S 1 whose 
generator is defined on cylinder functions by 

Gf(v) = E / e(^)ci(^»/)[/(^ ) -/(»/)]. (2.2) 

where / G = {/ : |||/||| = E ie z" A /« < °°} with A /« = sup{l/(»/) " /(C)l : 

^0) = CO') for a11 3 ^ *}• 

By Theorem 3.9 of Chapter 1 of Liggett (1985) the following condition, together with 
(|2.ip . implies that Q is the generator of a Feller process (o"t) on S : 

sup sup sup{ / p(da)\ci(a,rj) — Ci(a,r]^' b )\} < 00. (2-3) 

i<3L d »? fceA 7A 

In the following we shall work under conditions stronger than (|2 . 3[) ensuring not only that 
£/ is the generator of a unique Feller process, but also the possibility of perfectly simulating 
the stationary process corresponding to this infinitesimal generator. As a byproduct this 
implies that the system admits the existence of a unique invariant measure \x. 
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The main result of this article is a Kalikow-type convex decomposition of the change rates. 
We will prove that the change rate can be decomposed as 



Ci(a,rj) = Mi 



AiC-lJp^CoJ + ^AiC^Cal^WW)) 

k>0 



(2.4) 



where 

• Mi,i £ 1* d are positive constants, 

• for each i £ Z d , {Xi(k),k > -1} is a probability distribution on {—1,0, 1,2,...}, 

• for each i E Z d , p\ is a probability density on A with respect to the reference 
measure g, which does not depend on the configuration, 

• for each k > and for each r/ €E S, pj fe '(-|7/(Vi(fc))) i s a probability density with respect 
to the reference measure g, depending only on the local configuration r](Vi(k)). 



For convenience of the presentation we will add additional invisible jumps in (12. 2p . This 
is obtained by adding a cemetery A to A and defining A* := A U {A}. Define also 

g* := g + 5a- 

Denote 

Mi := sup / Ci(a,i])g(da). (2-5) 



Notice that M, is finite under condition (|2.ip . and define 

d{A, rj) := Mi - / Ci (a, r?)^(da). (2.6) 
J A 



mfa{A,rj) = 0. (2.7) 



Observe that 



Therefore we can rewrite the generator given by ()2.2[) as 

s/fa) = E / ^(^(".^[/c^ ) -/(»/)]. ( 2 - 8 ) 

where, by convention, for any z G Z d and any rj £ S = A^ d we define 

t/' A = 77. 

It follows that (|2.8p is a representation of the same generator as (12. 2|) . 

In order to obtain the decomposition we need the following continuity condition. 

Continuity condition. 

sup / sup \ci(a, rj) — q(o, Q\g{da) — > , (2-9) 

i&d J A v,CMVi(k))=C(V t (k)) 
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as k — > oo. 

To describe the convex decomposition of the rate function q, we have to introduce the 
following quantities. Define 

Qi(-l) = / inf Ci(a,C)g*(da), (2.10) 
J A* ceA ld 



and for any k > 0, 



aJk) = inf ( [ inf cj(a, C)e* (da) ) . (2.11) 



The continuity of Ci(a,rj) in 77 and the separability of 5 imply the measurability of 
inff:C(V i (fc))= u ,Ci(a,C) and inf CgAZd Cj(o,C) with respect to a. 

Note that by (j2?7j) 

/ inf Ci(a,()Q*(da) = / inf Cj(a, ()g(da). 
J A* c&A^ d ' .M ceA zd 

Further, by construction, we have that aj(fe) < on(k + 1), for each k > —1, and 

Mi = lim ai(fe). (2.12) 

k— >oo 

To obtain equality (j2. 12|) . fix some w € A Vi ( k > and observe that 
/ inf Ci(a,()Q*(da) = / inf Cj(a, Q^ida) 



+ Mj - sup / Ci(a,()Q(da). 
C4(Vi(k))=wJA 

/ inf Ci(a,()Q(da) — sup / Ci(a, ^)g(da) — > 

JAC-C(Vi(k))=w f.fty (k))=w J A 



as k — y 00 thanks to condition (|2.9p . 

Hence to each site i we can associate a probability distribution Ai by 



Ai(-l) = 4t^> ( 2 - 13 ) 



Mi 



and for k > 



Mfc) = Q ' (fc) ~ a ' (fc " 1) . (2.14) 
n ; Mi v ' 



Now we are ready to state the decomposition theorem. 

Theorem 1 Let (ci) i£ ^d be a family of measurable rate functions satisfying conditions 
\2.1\) . 112. 6)) and \2. 9\) . Then, for each site i, for Mi defined by 12. 5\) . and Aj(-) defined by 
\2. 13\) and {2.11$ , there exist 

• p\ a probability density with respect to g with support A, 
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• a family of conditional probability densities p\ (given by \3.21\) ), k > on A*, with 
respect to g* , depending on the local configurations r](Vi(k)) £ A Vi ^ 



such that 



Ci(a,rj) = MiPi(a\rj) for g* — almost all a € A* 



(2.15) 



where 



pMv)=H-m 1 \a) + Y,Hk)pf\a\7 ] (V l (k))). 



(2.16) 



As a consequence, the infinitesimal generator Q given by \2. 8\) can be rewritten as 



Note that for k = — 1, p\ (a) does not depend on the configuration and Aj(— 1) represents 
the spontaneous self-coloring rate of site i in the process. We will see in the proof that 
p\ 1 ' is defined in such way that p\ ^(A) = and therefore, the choice k = — 1 implies 
always a choice of a real color a € A, not of a = A. 

The decomposition given in Theorem [T] was designed in such way that the probability of 
self-coloring is maximized. This is important to speed up the perfect simulation algorithm. 
Obviously, slight modifications can be employed for different purposes as we will see in 
Example 2 (Section Hj). 

The representation given by (|2.17p provides a random finite range description of the time 
evolution of the process. We start with an initial configuration rj at time zero. For 
each site i E we consider a rate Mj Poisson point process N" 1 . The Poisson processes 
corresponding to distinct sites are all independent. If at time t, the Poisson clock associated 
to site i rings, we choose a range k with probability A,(fc) independently of everything else. 
And then, we update the value of the configuration at this site by choosing a symbol a 
with probability p^\a\at(Vi(k))) g* (da) . Choosing the symbol A means that we actually 
keep the current value of the spin. 

In Section we give examples of infinite range interacting systems where Theorem [T] can 
be applied. 




(2.17) 




fc>0 



3 Proof of Theorem [T] 



Put for any a £ A* 




inf Ci(a,(), 




c p(a\ v (i)) 

Af\a\n(i)) 





c^(a\n(i))-c^{a). 
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For any k > 1, define 

cf ] (a\ V (V t (k)))= inf a(a,C), 

AM(a|»7M(fc))) = cfVtoCW))) - cf^a^^ - 1)))- 
Then we have that for any a £ A, 

k 

Cl (a, V ) = *h*\vm))) + [hM - 4%\v(Vi(k))j\ . (3.18) 
i i 

Note that 

c]~ 1] (A) = inf Ci(A,r]) = Mi — sup / Ci(a,rj) 6 (da) = 0. 

Therefore, for a = A decomposition (|3.18p starts with j = 0, 

fe 

Ci(A,T/) = ^Al j] (A|ry(K(j))) + [ Q (A, ry) - <f(AM^(fc))) . 
i=o 

By monotonicity, we have for £* — almost all a € A* that 

cf 1 (a|7?(F i (/c))) -> limc| fe] (a|?7(V;(fc))) < q(o,t/) as fc -> oo. 
k 

Hence by monotone convergence, 

/ c? ] (a\r,(Vi(k))) e (da) -> / [limcf ] (a|ry(^(fc)))]^a) < / Cj (a,r?)^a). 

On the other hand, by (p^|) . 

cf ] (a|7?(yj(fe)))£>(do) -4 / Ci(a,rj)g(da). 
A J A 

Hence for g*— almost all a and all rj, 

oo 

limcf ] (a|rX^))) = ^ ] (a\r 1 (V i (j)))=c i (a, V ). 

Taking into account (|XTU|) and ([2TT5]) 

MiAi(-l)= / A!~ 1] (a)0(da). 

Hence we can define 



A 



Pi (a) " 



MiAi(-l) 
and 

p1- 1] (a) = o. 



Hence, p\ (a) is a probability density with respect to g*. Now, for k > 0, put 

Ai(M(W)) = ^ ^ Af^al^^)))^^), (3.19) 
and for any i, k such that Xi(k, T)(Vi(k))) > 0, we define 

Mi A i (/c,r ? (^(A:))) 

For i,A; such that Aj(/c, r/(Vi(fc))) = 0, define p^i\a\r]{Vi (k))) in an arbitrary fixed way. 
Hence for g* — almost all a £ A*, 



Ci(a,7]) = Mi 



\ i (-l)pl 1 \a) + J2HkMV i m)P?\a\v(Vi(k))) 



fc=0 



(3.20) 



In (|3,20p the factors Xi(k,r](Vi(k))),k > 0, still depend on r/(Vi(k)). To obtain the decom- 
position as in the theorem, we must rewrite it as follows. 

For any i, take Mi as in (|2.12p and the sequences a«(fe), Xi(k),k > —1, as defined in (|2. 1 1 j) 
and (|2.14p . respectively. Define the new quantities 

a i (kMVi(k)))=M i J2~ x i(l,r j {V i (l))) 

l<k 

and notice that 



O:, 



{k^Viik))) = / ^ ] (a,r,(Vi(k)))g*(da) 



A* 



is the total mass associated to cf : \. : r](Vi(k))). 

By definition of ai(k) in (|2.1ip . ai(k) is the smallest total mass associated to uniformly 
with respect to all possible neighborhoods rj(Vi(k)). Hence, in order to get a decomposition 
with weights Xi(k) not depending on the configuration, we have to define a partition of 
the interval [0, ai(k, rj(Vi(k)))] according to the values of ai(k). 

This yields, for any k > 0, the following definition of the conditional finite range probability 
densities. 

pf\a\ V (V i (k))) = 
fc-i 

^2 l {a l {l'-l,r,(V l (V-l)))<c il {k-l)<aS>,r l {V l ^ 
-l=l'<l 

^ ^ MiXi(k) P * 

m=l'+l 1 n ' 



+ a ' <fc> ^titr ( " )) ^ + " ( '"" w( ' +i))) 

The desired decomposition now follows from this. 



(3.21) 
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4 Examples 



In this section we give examples where the decomposition of Theorem [TJ can be applied. 
We start with an example from Bayesian statistics and image reconstruction. 



Example 1. Autonormal distribution. This model can be seen the spatial analogue 
of the autoregressive model. The usual way to describe its dynamics is through the si- 
multaneous schemes: Each pixel updates its value using a normal distribution with mean 
depending on the values of its neighbors. 

In this work, we are going to generalize this definition to incorporate long range interac- 
tions and arbitrary neighborhoods, but we are going to limit the values of the process to 
be on a compact interval [l,u]. The existence of a unique invariant measure for this pro- 
cess was studied by McBryan and Spencer (1977) and revisited by Ferrari and Grynberg 
(2008). Gibbs (2004) proposes a finite version of this model as the posterior distribution 
for Bayesian restoration of greyscale images. Huber (2007) studies perfect simulation of 
these distributions in a finite box. 

Consider that each pixel has an independent exponential clock and when the clock rings 
at the pixel i it updates its value depending on the values of its neighbors using a normal 
distribution with mean h(i,rj) and variance cr 2 (i) restricted that its value belong to a 
compact interval, without loss of generality we can consider the interval [0, 1]. The term 
h(i, rj) depends on r\ only through the values on the neighborhood of the site i. Typically, 

where J : % d — > R + is summable, non-negative and symmetric: J(i) > for all i € 
J(0) = and < J := ^2 i&d J(i) = 1. 

In this case, g(da) = \o,i] (a) da and 

/ \ 1 0((a - h(i,r]))/a(i)) 

a(a,ri) = — — — — — — — — ; — — — — -— , < a < 1 4.22 

where <ft and are the density and cumulative function of the standard normal distribution, 
respectively. 

Applying the bounds given by Proposition 1 of Fernandez, Ferrari and Grynberg (2007) 
we can show that (|4.22p satisfies the assumptions (|2,ip and (|2,9p needed in order to apply 
the decomposition of Theorem [TJ In our case, Mj = 1 and 

«M,^(,(«^)-,(^)) + 

where fjT = inf^ h(i, rj), [i + = sup^ h(i, rj), 



a{i) J \ a(i) J ' \ a(i) J ly / ^W 



9 



and 

Xi{a{t)) = t - lo S-7+- 

2 fi + — fi 



Also, 



n 1 veAW Af(k)\ V *(«) / V *(*) yy 



Ar(k) V V y V ^(i) 

where ^~(/c) = inf^y^))^ h(i,rj), /if(k) = sup v , v{Vi{k))=w h(i,rj), 

at w = * f 1 -?:<*> 'j fz«sw i , Am = * f'-ffl 1 ) -• 1 - (fc> ' 



all) / \ an 



and 

*'(*."(•)) = — ^- lo g^y 

As second example we show that the decomposition presented in Theorem Q] can be effec- 
tively implemented in Gibbsian systems with compact- valued spins. We take A = [—1,1] 
and introduce the following definitions. 

Definition 1 A pairwise potential is a collection {J(i,j), (i,j) £ 7* d x Z d } of real numbers 
which satisfies 

J(i,i) = 0, sup \ J(hj)\ < °°- ( 4 - 2 5) 

i£Zd i6Z«« 

In what follows we use the notation 



For any i € Z d , let ry(i) be the value of the spin at site i in the configuration rj € 5. 

Definition 2 ^4 probability measure /i on (S,S) is said to be a Gibbs state relative to the 
potential {J(i,j)} if for all i € Z d , a version of the conditional probability density ofr](i), 
given r](j),j ^ i, is given by 



exp U£#i J(i,jMj) 
fj,(r}(i) = a\rj(j) for all j / i}) = 



where 



Z 11 = exp aVj(g>(i) g(da). 
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In the following we consider the interaction Jr = (3J, where j3 is a positive parameter. The 
associated Gibbs measure will be denoted without indicating explicitly the dependence 
on j3. Now, put 

c i (a,7 ] ) = e^ a ^ dJ{i ' jUj) . (4.26) 

Then, by construction, the process (at) with generator (|2.2[) and this choice of change rates 
is reversible with respect to the Gibbs state fi corresponding to the potential Jp(i, j) = 
/3J(i,j). It is immediate to see that condition (|4.25|> implies the continuity condition (12. 9p . 

We now give the explicit decomposition in one specific case. 



Example 2. The following example is a Gibbsian time evolution with infinite range 
interaction. The decomposition we present here is inspired by the one presented in Galves, 
Locherbach and Orlandi (2010) in the case of two colors systems. In Galves et al. (2010), 
for coupling reasons, it was convenient to give a slightly different decomposition. The goal 
there was to be able to couple together the infinite range Gibbsian system with the finite 
range Gibbsian system obtained by truncating the potential interaction. We suppose that 
the spin distribution g is symmetric. 

Define for any i € Z rf and any k > — 1, 

S> k := \J(i,3)\, S^ k := £ \J(i,j)\. 

3 : \\i~ j : 1 1 i — j'||<fc 

>-l 



Note that £j = S % 
Then the decomposition (12.160 holds with 



M i= I (e a/3S * +e" a/3S M g{da). (4.27) 



Moreover, 



and 



Finally, 



and 



o 



o i (-l) = 2/ e- a ^'g(da) (4.28) 
Jo 



ou(k) = Mi + 



f 1 e ^ k e -aPS>X g{da) _ f 1 e a^ e(da y (4 29) 

Jo Jo 



Ai -1 = 2—1 — *^_J 4.30 

f (e a ^+e- a ^)g(da) 



1 papS^- 1 e -apS. 



Jo* 



> k ' e a P E 3 -:||i-i||=fc \J(h3)\ 



-«0Ej:||i-i||= h |J(<J)l 



g(da) 



(4.31) 
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5 Perfect simulation 



The goal of this section is to give an application of the Kalikow-type decomposition given 
by TheoremQ] This application is a perfect simulation algorithm for the invariant measure 
of an interacting multicolor system. We assume that the interaction rates are continuous 
in the sense of (|2,9p and satisfy a high noise condition. The basis of the algorithm is the 
convex decomposition given in Theorem [TJ First of all, the Proposition [1] gives a sufficient 
condition for exponential ergodicity which is based on the construction of a dominating 
branching process. 

From now on we will denote by (er^) (and (c^)) the multicolor system having generator 
Q given by (|2,8p with a fixed initial configuration n (a random configuration chosen with 
probability distribution /j,). 

Proposition 1 Let (ci) ie %d be a family of rate functions satisfying the conditions of The- 
oremUl Furthermore, assume that 

M= inf Mi > (5.32) 

and 

sup^ \Vi(k)\Xi(k) = 7 < 1. (5.33) 

Then the following two statements hold. 

1. The process (at) admits a unique invariant probability measure /i. 

2. For any finite set of sites F C Z d , for any T > and any initial configuration rj, 
there exists a coupling between the process (a^) and the stationary process (a^) such 
that 

P(4(F) + a£(F)) < \F\e-^-^ T . 

Let us compare the above proposition to known results in the literature on particle systems. 

1. Condition f|5.33|) is stronger than Liggett's existence condition (|2.3p which does not 
imply the uniqueness of the invariant measure. 

2. Let us compare our result to the M < e— criterion of Theorem 4.1 of Liggett (2000), 
page 31. Recall that Liggett's quantity M - translated into our context - is given 
by 

M = sup SU P {II c i (a> rj)g(da) - Ci(a,()e(da)\\ T v ■ r](k) = ((h) Vk / j} . 

By our decomposition of Ci(a,n), this expression can be upper bounded by 

M< sup Mi J2HkM(k)\. 

Since sup. Mi < oo, condition sup igZ d ^ fc>0 \Vi(k)\Xi(k) < oo implies Condition 
(3.8) of Liggett (2000). 
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Concerning the quantity e defined on page 24 of Liggett (2000), note that in our 
case, it can be written as follows. 



e= inf inf \ Cl (a,^ b )g({a}) + Cl (b, V l ' a )g({b})], 

which will be equal to zero in general. Hence we are able to treat with our techniques 
cases where the M < e— condition of Liggett (2000), Theorem 4.1, is not satisfied. 

3. Condition (|5,33p is a high-noise condition reminiscent of Dobrushin-Shlosman condi- 
tion, see Maes and Shlosman (1991). It is a sufficient condition ensuring that there 
is no phase transition. 

We are now in position to present the perfect simulation scheme. Suppose we want to 
sample the configuration at site i under the invariant measure \x. In a first step, we deter- 
mine the set of sites whose spins influence the spin at site i under equilibrium. We call 
this set of sites ancestors of % and this stage backward sketch procedure. First, we climb up 
from time using a reverse time Poisson point process with rate Mj. We stop when the 
last Poisson clock before time rings. At that time, we choose a range k with probability 
Xi(k). If k = — 1, we decide the value of the spin using the law p\ ^dg, independently of 
everything else. If k is different from —1, we restart the above procedure from every site 
j G Vi{k). The procedure stops once each site involved has chosen range —1. When this 
occurs, we can start the second stage, in which we go back to the future assigning spins to 
all sites visited during the first stage. We call this procedure forward spin assignment pro- 
cedure. This is done from the past to the future by using the update probability densities 
p\ starting at the sites which ended the first procedure by choosing range —1. For each 
one of these sites a spin is chosen according to p^~^dg. The values obtained in this way 
enter successively in the choice of the values of the spins depending on a neighborhood of 
range greater or equal to 0. 

We now give the precise form of the algorithm. Fix a finite set F C The following 
variables will be used. 

• N is an auxiliary variable taking values in the set of non- negative integers {0, 1,2,...} 

• Ngrp'op is a counter taking values in the set of non-negative integers {0, 1, 2, . . .} 

• I is a variable taking values in 7j d 

• K is a variable taking values in {— 1, 0, 1, . . .} 

• B is an array of elements of Z rf x {—1, 0, 1, . . .} 

• C is a variable taking values in the set of finite subsets of Z d 

• W is an auxiliary variable taking values in A* 

• a is a function from Z d to A* 

Algorithm 1 Backward sketch procedure 
1. Input: F; Output: N^Qp, B 
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2. N <- 0, N^op <~ °> B <~ 0> C <- ^ 

3. WHILE C / 

4. iV <- A + 1 

5. Choose randomly a position I E C and an integer K > — 1 according to the proba- 
bility distribution 

6. IF K= -1, C <- C\ {/} 

7. ELSE C^-CU Bj(K) 

8. ENDIF 

9. B(iV) «- (J, if) 

10. ENDWHILE 

11. N$ Qp <- iV 

12. RETURN N s ^ op , B. 

Now we use the following Forward spin assignment procedure to sample from the invariant 
measure [i. Recall that the choice of A in (|2.8p implies that the system does not change 
its colors. This explains Step 9 in Algorithm 2. 

Algorithm 2 Forward spin assignment procedure 
1. Input: Ngp p, B; Output: {(i,a(i)) : i G F} 

Z. 1\ <r~ JV STOP 

3. <- A for all j € Z d 

4. WHILE N >1 

5. (J, if) <- 5(iV). 

6. IF K = — 1 choose randomly in ^4 according to the probability distribution 

p\ l ^dg 

7. ELSE choose W randomly in A* according to the probability distribution 

pf\-\o)d B * 

8. ENDIF 

9. IF W ^ A put a(I) <- W 
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10. ENDIF 



11. N <- N - 1 

12. ENDWHILE 

13. RETURN {(«,<r(i)) :«ef) 

The next theorems summarize the properties of Algorithms 1 and 2. In order to distin- 
guish clearly to which part of the two algorithms we refer, we shall write P sketch for the 
probability associated to the backward sketch procedure. 

Theorem 2 Suppose that 

for all FcZ d finite, P a ketch(N$ OP < oo) = 1. (5.34) 
Then the following two statements hold. 

1. Algorithms 1 and 2 are successful. 

2. The process (at) admits a unique invariant measure [i. The law of the set {(£, : 
i £ F} printed at the end of Algorithm 1 and 2 is the projection on A F of \x. 

The next theorem states sufficient conditions ensuring (|5.34p and gives also a control on 
the rate of convergence. 

Theorem 3 1. The sub-criticality condition \5. 33\) implies {5.3$ - More precisely, we 
have 

Psketch(N { s % P >N)< \F\j N , (5.35) 

where 7 is given in \5.33\) . 

2. Suppose in addition to \5. 33\) that \5.32}) holds. Fix a time t > 0, some finite set 
of sites F C 1j d and two initial configurations n and £ € A^ d . Then there exists a 
coupling of the two processes (<Js)s an d (<Js)s such that 

P(a?(F) * ^(F)) < |F|e-M(i-7)*. 
The proofs of Proposition [lj Theorem [2] and Theorem [3] will be given in the next section. 

6 Proofs of Proposition CD, Theorems [2] and [3] 

The proofs rely on the notion of black and white time-reverse sketch process that we will 
introduce now. The black and white time-reverse sketch process gives the mathematically 
precise description of the backward black and white Algorithm 1 given above. 

For each i 6 Z d , denote by . . . T % _ 2 < Tl 1 < Tq < < T\ < T\ < . . . the occurrence times 
of the rate Mj Poisson point process on the real line. The Poisson point processes 
associated to different sites are independent. To each point T^ associate an independent 
mark K l n according to the probability distribution (Aj(A;))fc>_i. As usual, we identify the 
Poisson point processes and the associated counting measures. 
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For each i £ Z rf and t G R we define the time-reverse point process starting at time i, 
associated to site i, 

^(.ht) f _ rpi + > n 

1 n — L 1 iv*(0,t]-n+l> 1 — u ' 

= t-r Ni(m _ n+1 , t<0. (6.36) 

To these time-reverse point processes we can associate in an obvious way the corresponding 
marks Kn ,l \n € Z. Finally, for each site i € Z^, k > —1, the reversed fc-marked Poisson 
point process returning from time t is defined as 

iV ( ^[^]=El {# <, } l { ^) =fcr (6-37) 

To define the black and white time-reverse sketch process we need to introduce a family of 
transformations {TT^ ,k \i 6 Z d ,k > —1} on the set of finite subsets of Z d , F(7j d ), defined 
as follows. For any unitary set {j}, 

= {Sf 

Notice that for /c = —1, 7r^ ,fc - ) ({z}) = 0. For any set finite set F C Z d , we define similarly 

tt^(F) = U i6J r7T^({j})- (6-39) 

The black and white time-reverse sketch process starting at site i at time t will be denoted 
by (C'i*'*' > ) s >o. Cj ' is the set of sites at time s whose colors affect the color of site i at 

(i t) 

time t. We call this set Cg ' set of ancestors of i at time s before time t. The evolution of 
this process is defined through the following equation: Cq := {i}, and 

/(Cf>*)) = /(C£ M) ) + E E r[/(^' ,fe) ri- ) ))-/(^ ) )]^ (i ' t ' fc) (^) ; (6-40) 

where / : T(7L d *) — >■ M is any bounded cylindrical function. This family of equations 
characterizes completely the time evolution {ci l ''\s > 0}. For any finite set F C Z d 
define 

The following proposition summarizes the properties of the family of processes defined 
above. 

Proposition 2 For any finite set F C Z d , {Cg *\ s > 0} is a Markov jump process 
having as infinitesimal generator 

Lf{C) = M 8 ^ ^ A,(fc)[/(C U F 8 - /(C)] + ^(-^[/(C \ - /(C)], (6.41) 

i€C k>0 

where f is any bounded cylindrical function. 

Proof The proof follows in a standard way from the construction (|6.40p . • 

If we are interested in simulating from the invariant measure of the process, then we will 
start the black and white time-reverse sketch process at time t = 0, if however we wish 
to construct the process at time t, we shall start the black and white time-reverse sketch 
process at that time t precisely. 
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6.1 Backwards oriented percolation and sub-criticality 



For the algorithm to be successful it is crucial to show that U s >o C« , the set of ancestors 
of any site i, is finite with probability one. Formally, let 

TH OP = inf{, : CP = 0} 

~ (i) 

be the relaxation time. We introduce the sequence of successive jump times T„ , n > 1, 
of processes N^' k ^ whose jumps occur in (^01) . for t = 0. Let ff 1 = lf' 0) and define 
successively for n > 2 

fjp = inf{t > T«! : 3j € eg? , 3fc : A^QT^, *]) = 1}. (6.42) 

We write iT„ f° r the associated marks. Now we put 

C« = eg? (6.43) 

and 

A# OP = inf{n:Cf =0}. 

This is the number of steps of the backward sketch process - and it is exactly the number 
of steps of Algorithm 1. For the perfect simulation algorithm to be successful, it is crucial 

(i) 

to show that the number of steps NgJp OP is finite. Since at every step of the algorithm 

(i) 

a finite interaction range k is chosen, this implies automatically that also Tgj, OP is finite 
almost surely. However, in order to control the speed of convergence, we need a precise 

(i) 

control on the tail probabilities of T^ 0p . To this aim we estimate the volume of the set 
C { /' t} = U ieF cf '* } where F is a bounded set of 7L d . 

Lemma 1 

E(\Cf'% < |F|e" M(1 - 7)s ' (6.44) 
where M_ is defined in (|5.32p and 7 in (|5.33[) . 

Proof Fix some JVeN. Let L\ = |cf '* } | and 

= inf{t : L\ > N}. 

Then by dOD]) . 

rsAT N 



l\ ATn < i+EE/ [mk)\-^ {3&c ^ } ^' m {du) 

k>ijez d 

rs/\T N 

-E L W^ -1 ^- (6 - 45) 

Recall that that M_ = inf igZ d Mi > 0. Passing to expectation and using that, by condition 
M, ( C£ ^(kWm\ - 1]) - A,(-l) J < -Mil - 7) < 0, 

V fc>i / 
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this yields 



e(l\ ATn ) < i+x; M i|(E A /W[iwi-i])-^(-i) 

fc>i 

l {jeC ,M )} dn 

< 1 - M(l - 7) E I V u du. (6.46) 



Letting N — > 00, we thus get by Fatou's lemma that 

£(4)<l-M(l- 7 ) I* E{I> u )du. 

Jo 

This implies that 

< 1 for all s > 0. 
Hence we may apply Gronwall's lemma which yields 

E(L\) < e -Mi-7)s. ( 6 .4 7 ) 

Hence, since \ci F,t) \ < ^ ieF \ c s\ = Y,ieF L l> 

E(\C^\) < iFle-^ 1 ^ 3 - (6.48) 



6.2 Proof of Proposition [T] and Theorem [2] 

Proposition [1] is an immediate consequence of Theorem item 2, and Theorem El item 2. 
Proof of Theorem [2j 

Item 1 of Theorem [2] is evident. We give the proof of item 2 of Theorem [2j 

Write fiF for the law of the output {(i,a(i)) : i £ F} of Algorithm 1 and 2, \ip is a 
probability measure on (A F , A F ). By construction, the family of probability laws {/ip, F C 
Z d finite} is a consistent family of finite dimensional distributions. Hence there exists a 
unique probability measure fi on (S,S) such that \ip is the projection onto A F of /i, for 
any fixed finite set of sites F C 'L d . 

We show that [i is the unique invariant measure of the process (at)- In order to do so, we 
use a slight modification of Algorithm 1 and 2 in order to construct crj 7 , for some fixed 
initial configuration 7/ G S. The modification is defined as follows. Let Tstop an d T be 
variables taking values in (0, 00). Replace Steps 1 — 3 of Algorithm 1 by 

1. Input: F] Output: N^ OP ,B,C. 

2. N <- 0, iVlrop <- 0, £ <- 0, C <- F, T S top <- 

3. WHILE Tstop < t and C ^ 
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3'. Choose a time T G (0, +00) randomly according to the exponential distribution with 
parameter Yljec^j- Update 

Tstop Tstop + T. 
Finally replace Step 12 of Algorithm 1 by 

12. RETURN N$ OP , B > C. 

In this modified version, we stop the algorithm after time t, hence the output set C might 

(Ft) 

not be empty. The output C is exactly the set C\ , the set of sites at time whose 
colors influence the colors of sites in F at time t. Finally, notice that if C = 0, then 
Tstop < t, and in this case, Tstop = T F TOp is the relaxation time introduced in the 
previous subsection. 

Concerning Algorithm 2, replace Step 1 of Algorithm 2 by 
1. Input: N$ op ,B,C; Output: {(i,<r(i)) ■ i G F}. 
and Step 3 by 

3. a{j) <- 7](j) for all j G C; a{j) <- A for all j€Z d \ C. 

Then the law of the set {(i, a(i)) : i G F} printed at the end of the modified Algorithm 2 
is the law of a^(F). Notice that the output of the modified Algorithm 2 equals the output 
of the unmodified Algorithm 2 if Tstop < t. 

We first give an intuitive argument showing that fi must be invariant for (er t ). Write Pt 
for the transition semigroup of (at). Fix t > and a finite set of sites F C U 1 . Suppose 
we want to determine the projection on A F of fiPt. This means that we have firstly to 
run the above modified Algorithm 1 up to time t. It gives as output the set of ancestor 
sites C = C\ . We then have to run the modified Algorithm 2 with initial configuration 

( F t) 

{a(i),i G C = C t } chosen according to fi (cf. Step 3.). But this means that we have to 
concatenate the modified Algorithm 1 with the original Algorithm 1, where the Algorithm 

1 is now starting from C^'*" 1 instead of F. In other words, we concatenate two backward 
sketch processes and consider ' > up to the time of its extinction. By the Markov 
property and the stationarity of the Poisson processes, this means that we consider directly 
C( F >tf up to the time of its extinction - which is finite by our assumptions. Hence \iPt = 11 
in restriction to finite cylinder sets. 

We now give a more formal argument. Recall that the output of the modified Algorithm 

2 equals the output of the unmodified Algorithm 2 if Tstop < t. Let / : A F — > R + be a 
bounded measurable function. Then 

E[f(a r l(i), i e F)] = E[f(a?(i),i G F),T STO P < t] + E[f(a?(i),i G F),T STO P > A 
= E[f(a(i),ieF),T F TOP <t] 

+E[f(a?(i),i G F),T S top > t], (6.49) 

where (a(i),i G F) is the output of the unmodified Algorithm 1 and 2. 



19 



But 

E[f(<r v t (i),i G F),T STOP >t}< \\f\\ooPsketch(Tg TOP > t) -> as t -)■ oo, 
since finiteness of N^ TOp implies the finiteness of Tg TOp . Hence we obtain that 

lim E[f(^(i),i G F)j = E[f(a(i),i G F)], 

i— >oo 

since 1{tJ to <t} ~~ ^ ahnost surely. 

This implies that jx is an invariant measure of the process. Replacing the initial condition rj 
by any stationary initial condition, we finally get also uniqueness of the invariant measure. 
Thus Theorem [2] is proved. • 



6.3 Proof of Theorem d 

We start by proving f|5.35j) . Let 

L (i) = | C (i)| 

be the cardinal of the set c4 after n steps of the algorithm (recall (I6.43P ). Then due 

(i) 

to our assumptions, L n can be compared to a multi-type branching process Z n having 

(i) 

offspring mean which is bounded by 7 at each step, such that Ln < Z n for all n. The 
details are given in Galves, Locherbach and Orlandi (2010), proof of Theorem 1, and are 
omitted here. Thus, 

P(N^ Qp >n) = P(L« > 0) = P(LW > 1) < P(Z n > 1) < E{Z n ) = 7 ». 

When starting with the initial set P instead of the singleton {i}, then the above estimates 
remain true by multiplying with |P|, due to the independence properties of the branching 
process. 

Concerning item 2 of Theorem [3l we use once more the modified algorithms 1 and 2 
introduced in the proof of Theorem [2] above. In order to realize the coupling, we use the 
same realizations of T, I and K for the construction of and a^. Write L s for the cardinal 
of d F,t \ Clearly, both realizations of erj' and do not depend on the initial configuration 
rj, £ respectively if the output C of the modified Algorithm 1 is void. Thus, by Lemma [U 

P(aUF)^a C t (F)) < P(T ST OP>t) 
= P(L t >l) 
< E(L t )<\F\e-^ 1 -^. 

This implies that the convergence towards the unique invariant measure takes place expo- 
nentially fast. This concludes the proof of Theorem [3J • 



7 Applications for perfect simulation 

7.1 Maximum likelihood estimation in Gibbs distributions 

Parameter estimation for Gibbs distributions in the infinite lattice is usually based on 
the maximum likelihood approach, see for example Gidas(1988, 1991). The maximum 
likelihood estimation is theoretically well understood in this framework. Comets (1992) 
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proved the consistency of the MLE for exponential families of Markov random fields on 
the lattice. Also, in the case of no phase transition, Janzura (1997) proved asymptotic 
normality and efficiency of the MLE inside the uniqueness region of the Gibbs distribu- 
tions considered. Comets and Gidas (1992) considered maximum likelihood estimators for 
Markov random fields over Z d from incomplete data. They prove the strong consistency 
of maximum likelihood estimators in this case. Their results hold irrespectively of the 
presence of long-range correlations or non-analytic behavior of the underlying quantities. 
The parameter space is thereby allowed to be non-compact. 

However, the numerical feasibility of the ML method is strongly limited, due to the com- 
putation of the normalizing constant for each relevant parameter, in particular for each 
temperature. Geyer and Thompson (1992) devised a rather ingenious method for this 
computation based on an MCMC computation of the equilibrium distribution for a fixed 
value of the parameter. 

Consider the family of probability densities with respect to a reference measure \x given 
by 

f(x,0) = ^exp<T(x),6> (7.50) 

where < T(x),9 > denotes the inner product between the canonical parameter 9 and the 
sufficient statistics T{x) and 

Z(6) = J exp < T(x),9 > dfx(x). (7.51) 

Denote by P^ the measure having density f(-,ip) with respect to fi. Then, 

Z{9) = Z(i/>) J exp < T(x),6- ip > dP^(x). (7.52) 

Therefore if we have X±, X2, ... iid random objects with distribution P^ we have that 

1 n Z(6) 
d n {9) = - exp < T(Xi), 9-ip>^ d{6) = almost surely. (7.53) 

The maximum likelihood of 9 can be taken as 

= arg maxlog f(x, 9) + log Z(tp) = arg max < T(x), 9 > — log d(9), (7.54) 
and its Monte-Carlo approximant 

9 n = arg max < T(x),9 > - log d n (9). (7.55) 

Notice that if we can perfectly simulate from Pw, we have trivially that 9 n — >• 9 as n — > 00 
along with the rate of convergence of such convergence. 



7.2 Attaining Ornstein's d— distance for ordered pairs of Ising probabil- 
ity distributions 

In this section we show how to use the decomposition of Theorem [T] and the above per- 
fect simulation algorithm in order to construct an explicit coupling attaining Ornstein's 
J— distance for two ordered Ising probability measures. Let A := {—1, 1} and S = A z . 
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We consider a ferromagnetic pairwise interaction J, i.e. a collection {J(i,j),i ^ j,i,j G 
7L d } of positive real numbers satisfying J(i,j) = J(j,i) for all i,j G Z d and for all i G Z d , 

J(i,i)=0, sup J(i,j) < oo. (7.56) 

Let {J(i,j),i ^ j, i, j G Z d } be another pairwise interaction satisfying an analogous 
summability assumption such that 

J(i,j)<J(i,j)ioi all i,jeZ rf 

Moreover, let {h{,i G Z '} and G Z d } be two collections of positive real numbers 

representing an external field such that 

hi < hi for all i G Z d , sup hi < oo. 

i 

Finally we suppose that 

for all i G Z d , ^[J(i, j) - J(i,j)} <hi-h % . (7.57) 

3 

Recall that a probability measure \i on S is said to be a Gibbs measure relative to the 
interaction J and the external field {hi} if for all j G Z d and for any fixed £ G 5, a version 
of the conditional probability /x({cr : a(i) = ((i)\a(j) = C(j) f° r & U J 7^ *}) ^ s gi ven by 

fi({a : a(i) = ((i)\a(j) = ((j) for all j ^ »}) = — 1 rr-r— - tt^t . 

1 + exp(-2/3[X; j J(«,j)CWCu) + ^CWD 

(7.58) 

The Gibbs measure /2 associated to the interaction J and the external field {hi} is intro- 
duced analogously. 

We consider a Glauber dynamics (at(i),i G G R) taking values in S = A zd and 

having \i as reversible measure. The process is defined by the rates Ci(a), i G Z d , where 
Cj(cr) is the rate at which the spin i flips (i.e. changes its sign) when the system is in the 
configuration a. We take 

d(a) = exp J(hj)<?(i)v(j) + M*)] j ■ (7-59) 

By construction, the process (cr^ is reversible with respect to the Gibbs measure fx. In 
the same way we can define a Glauber dynamics (&t)t reversible with respect to the Gibbs 
measure /}, associated to the interaction J and the external field {hi}. 

The main idea of our approach is a coupled construction of the processes (at) and (at) 
which is order preserving. More precisely, let us write 

a < a if and only if a(i) < a(i) for all i G if 1 . 

Write 

S = {(a, a) G A zd xA zd :a<a} = {(-1, -1), (-1, +1), (+1, +l)f d . 



22 



We now describe the coupled time evolution of a t and a t . Start with an ordered couple of 
initial configurations n < fj at time 0. Let 

M . _ 2e l 3 ^jezd J(.h3)+hi] m 

For each site i £ Z d , consider a rate Mj Poisson process iV\ The Poisson processes corre- 
sponding to distinct sites are independent. If at time t, the Poisson clock ate site i rings, 
then both processes try simultaneously to update their spin at site i. Process a replaces 
spin a(i) by — er(i) with probability 

Cj((Tt) 

Mi ' 

and the process a replaces spin a(i) by — a{i) with probability 

Mi ■ 

Hence we can introduce the following probability measures. For any configuration a with 
a(i) = —1 we put 

= ^-,Pi(-l\a) = l- Pi (+l\a). 
In the same way, for any configuration a with a(i) = +1, we put 

p»(-ik) = i^>^(+ 1 l CJ ) = i -w(-ik)- 

The same definitions hold for £j with obvious modifications. Then we have by construction 
and thanks to condition (I7.57P that 

Pi(+l\a) < pj(+l|cr), whenever a < a. 

This stochastic order makes it possible to construct a coupled Glauber dynamics (crt,&t)t 
taking values in the space of ordered configurations S. 

At each jump time t of one of the Poisson processes iV*, the ordered configuration (at, Of ) 
is replaced at site i by the ordered pair 

(+1,+1) with probability P i ((+l,+l)|(<7 t ,or t )) = ^(+1^), 

(-1,-1) with probability ^((-1, -l)|(f7 t , a t )) =p i (-l\a t ), (7.60) 
(-1.+1) with probability Pi((-1, +l)|(cr t , a t )) = Pi (-l\a t ) - Pi{-l\a t ) . 

Now, it is straightforward to show that under our summability condition (|7.56p on the 
interaction J and J and due to the boundedness of the force of the external field, the 
transition probability Pj satisfies the continuity assumption ()2.9|) . Hence, the decomposi- 
tion of Theorem [T] can be applied and yields the following (compare also to Theorem 3.3 
of Galves, Garcia and Prieur (2010)). 

Corollary 1 There exists a sequence of transition probabilities Pk,k > —1, such that for 
any pair of symbols (a, b) € {(—1, — 1), (— 1, 1), (1, 1)} and any ordered pair of configura- 
tions (a, cr)e5: 

oo 

P((a,6)|(a,a)= £ A#)pJ fe] ((a,6)|(a,cO(^))). 
fc=-i 
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As in Galves, Locherbach and Orlandi (2010), it can be shown that for this decomposition, 
a sufficient condition for (|5.33j) is 



BUpJ^(fc)| ]T J(i,j) ] <OC 

k \j-\\j-i\\=k 



and /3 < /3 C , where f3 c is solution of 

fc>l \ j-\\j-i\\=k J 

In this case, we can extend the ideas of Galves, Garcia and Prieur (2010) from the case of 
chains of infinite order to infinite range Gibbs measures: 

Our perfect simulation algorithm simulates two ordered configurations belonging to S 
according to the invariant distribution of the coupled Glauber dynamics. This yields an 
explicit coupling of the two Gibbs measures /i and ft,. Since this coupling is ordered, the 
very nice argument of the proof of Theorem 3.6 in Galves, Garcia and Prieur (2010) tells 
us that this coupling attains necessarily the d— distance 

d(fi, /}) = inf sup{P(cr(i) ^ a(i) : (a, a) is a coupling of /i and Jl }. 

i 

Hence our perfect simulation algorithm enables us to construct explicit couplings achieving 
this distance, and as far as we know the problem of finding explicit solutions was addressed 
only for finite volume Gibbs measures up to now. 



8 Impatient user bias 

Perfect simulation procedures, very often cannot be run until the algorithm stops either 
by limitations of time or limitations of buffer. In this section we give upper bounds for 
the probability of these two types of errors. 

According to our construction, the perfect simulation algorithm of /i presented in this 
article is a function F : [0, l] NxZ to S such that, if (U n ) n = (U n (i),i € Z d ) n is a sequence 
of i.i.d. families, indexed by Z d , of uniform in [0, 1] random variables, then for any site 

, («) (i) 

i (zTj , there exists a stopping time Ngj^ p, such that F depends only on the first Ngj^op 
families of (U n ) n , i.e. for any measurable B G A, 

L "STOP J 

(i) 

Note that Ngrp 0p is not the number of uniform random variables that have to be simulated 

(i) 

in order to sample from /i, this number will in general be considerably larger. Ngj, OP is 
the number of steps of the backward sketch procedure. 

A first kind of "impatient user bias" occurs whenever the user, for reasons independent of 
the algorithm, has to stop the algorithm after, say N steps maximal. In this case, we do 
not sample from /i, but instead sample from 



P 



FdU^j)),, (£/> e b \ n sIop < N 

"stop 
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By Proposition 6.2 of Fill (1998), compare also to section 6 of Ferrari, Fernandez and 
Garcia (2002), the error made above can be bounded by 

p(n^ op >n) ^ yv 



< 



1-P(N^ op >N) 1-7*' 
see Theorem [3] above. 

At each step of the backward sketch procedure, a range of order k is chosen, where k 
is, in general, not bounded from above. In practical situations, however, a user will be 
limited in the choice of the interaction range and will restrict the simulation to the choice 
of ranges bounded by a a certain upper bound L that he decided to fix in advance. More 
precisely, writing 

if := inf{Z f : K® > L}, 



the use will therefore sample from the measure 



P 



Fmm, ■ • • , (u ( i} (j))^ e B\{Nf TOP < n} n {if > r> TOP } 

ly STOP 



In order to control the error made induced by this "space-time impatient user bias", we 
have to control 

pyf < t% op ). 

Using arguments similar to Lemma 2 of Galves, LLocherbach and Orlandi (2010), this can 
be bounded by 
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9 Final comments and bibliographical discussion 

In this work we study the equilibrium measure of systems with infinite range interac- 
tions satisfying fast decay of the long range influence on the change rate and a certain 
subcriticality-criterion. For Gibbs random fields, this regime has traditionally been stud- 
ied via cluster-expansion methods which either rely on sophisticated combinatorial es- 
timations (Malyshev (1980), Seiler (1982), Brydges (1984)) or inductive hypothesis and 
complex analysis (Kotecky and Preiss (1986), Dobrushin (1996, 1996a)). 

This is not the approach we follow here. Our approach is probabilistic, based on an 
explicit construction of the dynamics and gives probabilistic insight into the structure of 
the stationary law of the process, without combinatorial or complex-analysis techniques. 
Let us stress that our approach is not an alternative to cluster expansions. It has a different 
regime of validity and different aims. 

Our construction is reminiscent of Harris' graphical representation for particle systems 
and it is similar in spirit to procedures adopted in Bertein and Galves (1977), Ferrari 
(1990), Van den Berg and Steif (1999), Ferrari, Fernandez and Garcia (2002), and Garcia 
and Marie (2006) among others. However, all these papers only consider particular mod- 
els, satisfying restrictive assumptions which are not assumed in the present paper. Our 
approach works for any infinite range continuous interaction under the only assumption 
of high-noise. 
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There are several techniques for perfect simulation of Markov processes. Among the most 
popular ones figures Coupling from the Past (CFTP) originally proposed by Propp and 
Wilson (1996) and applied to several special cases in a vast literature. A good review 
can be found in Kendall (2005). This kind of technique applies to invariant measures of 
Markov processes with finite coalescence time. One main point of the CFTP technique 
is that one has to be able to control the coalescence times uniformly with respect to all 
possible starting points. This is an issue that becomes particularly difficult in the case of 
"big" state spaces. The problem of large state spaces can be overcome for processes with 
certain monotonicity properties or for some specific cases. For example, for spatial point 
processes there is a vast literature on the subject, we point the works of Kendall (1998), 
Kendall and Thonnes (1999), Kendall and M0ller (2000) among others. 

In our case, we sample directly from a time stationary realization of the process. There is 
no coalescence criterion, either between coupled realizations or between sandwiching pro- 
cesses. The scheme neither requires nor takes advantage of monotonicity properties. The 
scheme directly samples a finite window of the equilibrium measure in infinite-volume. 
In contrast, other CFTP algorithms, for example Kendall (1997, 1998), focus on finite 
windows with fixed boundary conditions, and the infinite-volume limit requires an ad- 
ditional process of "perfect simulation in space". We point out Ferrari (1990), van den 
Berg (1993) and van den Berg and Maes (1994) have also proposed construction schemes 
for (infinite-volume) Gibbs measures of spin systems that can be easily transcribed into 
perfect-simulation algorithms. 

For continuous state spaces in systems with finite number of components, Cai (2005) 
proposes a non-monotone CFTP but as he points out "the detailed construction of the 
non-monotone CFTP algorithm is problem specific". Fernandez, Ferrari and Grynberg 
(2007) construct perfect simulation for random distributions supported onad dimensional 
box, in particular, multivariate normal distributions restricted to a compact set. Connor 
and Kendall (2007) show that for a large class for positive recurrent Markov processes it is 
always possible to perform CFTP, although not always feasible. In general, for interacting 
particle systems with continuous state spaces, it seems to be out of reach to apply CFTP 
successfully. A recent paper by Huber (2007) succeeded in using CFTP for a very specific 
case of a continuous autonormal system restricted to a finite box. 

The notion of random Markov chains was introduced explicitly in Kalikow (1990) and 
Bramson and Kalikow (1993) and appeared implicitly in Ferrari, Maass, Martinez and 
Ney (2000) and Comets, Fernandez and Ferrari (2002). 
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